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Abstract: Weighted averaged finite difference methods for solving fractional 
diffusion equations are discussed and different formulae of the discretization of the 
Riemann-Liouville derivative are considered. The stability analysis of the different 
numerical schemes is carried out by means of a procedure close to the well-known 
von Neumann method of ordinary diffusion equations. The stability bounds are 
easily found and checked in some representative examples. 
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1. INTRODUCTION 

The number of scientific and engineering prob- 
lems involving fractional calculus is very large 
and growing. The applications range from con- 
trol theory to transport problems in fractal struc- 
tures, from relaxation phenomena in disordered 
media to anomalous reaction kinetics of subd- 
iffusive reagents. Recently, a fractional Fokker- 
Planck equation has been proposed to describe 
subdiffusive anomalous transport in the presence 
of an external field (Metzler et ai, 1999; Barkai et 
al, 2000; Metzler et at, 2000). For the force- free 
case, the equation becomes the fractional partial 
differential equation(Balakrishnan, 1985; Scalas 
et al, 2004; Metzler et at, 2000; Schneider et 
al, 1988) 

— t) ^ oDI^^—u{x, t) (1) 

where oD\~'^ is the fractional derivative de- 
fined by the Riemann-Liouville operator, K-^ is 
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the diffusion coefficient and 7 G (0, 1) is the 
anomalous diffusion exponent. There are many 
analytical techniques for dealing with these frac- 
tional equations. But, as also with ordinary (non- 
fractional) partial differential equations (PDEs), 
in many cases the initial and/or boundary condi- 
tions and/or the external force are such that the 
only reasonable option is to resort to numerical 
methods. However, these methods are not as well 
studied as their non-fractional counterparts. 

In this communication, some numerical meth- 
ods for solving fractional partial differential equa- 
tions, which are very close to the well-known 
weighted average (WA) methods of ordinary (non- 
fractional) PDEs, are considered. It will be shown 
that the stability of the fractional numerical 
schemes can be analyzed very easily and effi- 
ciently with a method close to von Neumman's 
(or Fourier's) method for non-fractional PDEs. 

2. FRACTIONAL DISCRETIZATION 
FORMULAE 

Two main steps will be considered to build nu- 
merical difference schemes for solving fractional 



PDE's. In the first step one discretizes the ordi- 
nary difFcrcntial operators d/dt, d^/dx^ as usual 
(Press et ai, 1992; Morton et ai, 1994). This 
will be done in Sec. 3. In the second step, one 
discretizes the Riemann-Liouville operator: 
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where \t/h] means the integer part of t/h. This 
formula is not unique because there are many 
different valid choices for cj^"-* (Lubich, 1986). 
Let ui{z,a) be the generating function of the 



coefficients i.e., 
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If the generating function is 

oj{z,a) = {l-zr (4) 

then we get the backward difference formula of or- 
der p = 1 (BDFl). This is also called the backward 
Euler formula of order 1 or, simply the Griinwald- 
Letnikov formula. These coefficients are uj^jf* = 
(— 1)'^(^) and can be evaluated recursively: 
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The generating function for the backward differ- 
ence formula of order p = 2 (BDF2) is 
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is the generating function for the backward differ- 
ence formula of order p = 3 (BDF3). Generating 
functions for higher-order BDF formulae can be 
found in (Lubich, 1986; Podlubny, 1999). Another 
type of discretization formula is that of Newton- 
Gregory of order p (NGp) (Lubich, 1986) whose 
generating function is 
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where the coefficients f2„ are defined by 
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average method, the diffusion equation is replaced 
by: 



(m+l) (m) , ^ ri 
Uj =Uj ' + XS 



(1 - X)S 

T{x,t), 



Am) 



(m+l) 

i-1 



4"' 



2»; 



m+l) 



,M 

u 



+ 

(m+l) 



+ 



(10) 

where A is the weight factor, T{x,t) the truncation 
error (Morton et al., 1994), and S = DAt/{Axf. 
Similarly, the fractional equation is replaced by 
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where Uj = qD^ u{xj,tm)- Inserting Eq. (2) 
into Eq. (11), neglecting the truncation error, and 
rearranging the terms, we finally get the fractional 

WA difference scheme 
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For A = 1 we recover the fractional explicit 
method discussed in (Yuste et al, 2003; Yuste 
et al, 2004). The method is impHcit for A 7^ 1. 
For A = one gets the (fractional) fully im- 
plicit method, and for A = 1/2 the (fractional) 
Crank-Nicholson method. Because the estimates 
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J of u{xj,tm) are made at the times mAt, 
m = 1,2,..., and because the evaluation of 

^u(xj,t) by moans of (2) requires knowing 

it is 



u{xj,t) at the times nh, n = 0,1,2, 
natural to choose h = At. In this case, 
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Before tackling Eq. (12) seriously, one must first 
know under which conditions, if any, the integra- 
tion algorithm is stable. 



3. FRACTIONAL DIFFERENCE SCHEMES 

The notation xj = jAx, tm = mAt and 
u{xj,tm) = Uj^^ ~ C/j'"^ will be used with 

[/j™^ being the numerical estimate of u{x, t) at 
the point {xj,tm,). In the non-fractional weighted 



4. STABILITY ANALYSIS 

The stability analysis of the integration difference 
scheme (12) will be carried out by means of the 
method used in (Yuste et al., 2003; Yuste et 
ai, 2004). Following the von Neumann ideas, one 
studies the stability of a single suhdiffusive mode 
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Fig. 1. Stability phase diagram for the weighted 
average methods versus the weight factor A. 
The hne corresponds to Eq. (15). The square 
corresponds to the case shown in Fig. 5 and 
the star corresponds to the cases shown in 
Figs. 6 and 7. 

of the form t/j"^ = Cme''?^'^"=. This mode will be 
stable as long as the (m stay bounded for m — > 
oo. Proceeding as (Yuste et al., 2003; Yuste et 
al., 2004) one finds that a fractional WA method 
is stable as long as 1/S > 1/5'x where 

l/5x =2(2A-l)w(-l,l-7). (15) 

This is the main result of the present work. From 
(15) one sees that a WA method is stable for 
any value of S' if A < 1/2 because the generating 
function lo{z, 1 — 7) for z = —1 is positive (5' is 
always positive: see Eq. (13)). For A > 1/2, S'x 
is positive, and there exist values of 1/S' smaller 
than 1/Sx, so that the fractional WA methods 
are unstable for these cases. Figure 1 shows the 
stability phase diagram for the WA methods. For 
A = 1 one recovers the stability bound 1/Sx = 
2w(— 1, 1 — 7) for the fractional explicit methods 
discussed in (Yuste et al, 2003; Yuste et al, 2004). 
The stability bounds of these explicit methods 
versus the anomalous diffusion exponent 7 for 
several Ricmman-Liouville derivative discretiza- 
tion formulae are shown in Fig. 2. Note that the 
stability region shrinks when the order p of the 
discretization formula increases. 

The fully implicit numerical method (A = 0) and 
the Crank-Nicholson method (A = 1/2) have been 
considered, together with their stability, by some 
authors (Sanz-Serna, 1988; Lopez-Marcos, 2002; 
Lubiche< al., 1996; McLeanei al., 1993; McLeanei 
al, 1996). However, their analysis is far more com- 
plex than that presented here. Also, to the best of 
my knowledge, neither generic implicit WA meth- 
ods (arbitrary A) nor the explicit method (A = 1) 
have been considered previously. In this vein, it 
is interesting to note that an explicit method for 
solving the fractional diffusion equation written 
in the Caputo form has recently been proposed by 
Ciesielski and Leszczynski (Ciesielski et al., 2003). 
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Fig. 2. Stability bound for the explicit method 
versus the anomalous exponent 7. The lines 
BDFl, BDF2, BDF3 and NG2 correspond to 
the theoretical stability bounds of the explicit 
method when the BDFl, BDF2, BDF3 and 
NG2 discretization formulae of the Riemann- 
Liouville derivative are used. The circles are 
numerical estimates of the BDFl stability 
bound (Yuste et al, 2003; Yuste et al, 2004). 
The squares correspond to the cases shown in 
Fig. 3 and the star is the case of Fig. 4. 

5. CHECK OF THE STABILITY BOUND 

The stability of the fractional difference schemes 
of Sec. 3 will be checked here by applying them 
to solve Eq. (1) with the initial condition u{x,t = 
0) = .t(1 — x) and the absorbing boundary con- 
ditions u{0,t) = u{l,t) = 0. The exact analytical 
solution of Eq. (1) is easily found by the method 
of separation of variables: 

8 °° 1 
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(16) 

where E^^ is the Mittag-Leffler function (Metzler 
et al, 2000; Mainardi et al, 2000). In Fig. 3 
we compare this exact solution with the results 
provided by the BDFl explicit integration scheme 
for anomalous diffusion exponents 7 = 0.5, 7 = 
0.75, and 7 = 1 for time t = 0.5 and = 1. The 
values of Ax used were Ax = 1/10, 1/20, and 1/50 
with S = 0.33, 0.4, and 0.5, respectively. These 
values of S are marked by squares in Fig. 2. They 
arc inside the stable region, which is confirmed 
by the well-behaved numerical solutions shown 
in Fig. 3. However, in Fig. 4 the BDFl-explicit 
numerical solution for S = 0.37 and 7 = 1/2 
shows a characteristic unstable behaviour. This 
is the expected behavior because 1/0.37 is smaller 
than l/Sx = 2cj(-l, 1 - 7) = 22-t for 7 = 1/2. 
This case is marked by a star in Fig. 2. 

Figures 5,6, and 7 show the numerical integration 
results for the WA implicit method with A = 0.8, 

7 = 1/2 and two values of S. For S = 0.55 one has 
1/S > 1/5'x = 1-2 X 2^/^ (this case corresponds 
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Fig. 3. Solution of the subdifFusion equation (1) 
with absorbing boundary conditions u{0, t) = 
u{l,t) = and initial condition u{x,0) = 
x{l — x). The symbols correspond to the 
BDFl-explicit numerical solution and the 
lines correspond to the exact analytical so- 
lution. The solution is shown for the time 
t = 0.5 for 7 = 0.5, S = 0.33, Aa; = 1/10 
(triangles), 7 = 0.75, S = 0.4, Ax = 1/20 
(squares) and j = 1, S = 0.5, Ax = 1/50 
(circles). These cases are marked by squares 
in Fig. 2 



s 




Fig. 4. BDFl-explicit numerical solution u{x,t) 
for the same problem as in Fig. 3 but with 
7 = 1/2, ^ I, S ^ 0.37, and Ax = 
1/20 after 150 time steps (squares) and 200 
time steps (circles). The lines are plotted as 
a visual aide. This case corresponds to the 
point marked by a star in Fig. 2. 

to the square in Fig. 1) and the WA method must 
be stable. This is confirmed in Fig. 5. However, 
l/S" < l/^x = 1.2 X 2^/2 for S = 0.7 so that the 
WA method must be unstable in this case, which 
is confirmed in Figs. 6 and 7 (this case is marked 
by a star in Fig. 1). 



6. FINAL REMARKS 

The truncation error T{x,t) of Eq. 11 can be 

estimated in the same way as for non-fractional 
equations (Morton et al, 1994). One gets: 



Fig. 5. Numerical solution u{x,t) obtained by 
means of the weighted average method with 
weight factor A = 0.8 for the same problem as 
in Fig. 3 but with 7 = 1/2, K.^ = 1, S = 0.55, 
and Ax =1/20 after 500 time steps (circles). 
The line is the exact analytical result. This 
case corresponds to the point marked by a 
square in Fig. 1. 
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Fig. 



6. Numerical solution u{x,t) for the same 
problem as in Fig. 3 obtained by means of the 
weighted average method with weight factor 
A = 0.8 for 7 = 1/2, = 1, S = 0.7, 
and Aa; = 1/20 after 50 time steps. The 
lines are plotted as a visual aide. This case 
corresponds to the point marked by a star in 
Fig. 1. 
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where T™ is the truncation error at Xj and tm + 
At/2: Tp = T{xj,tm+i/2)- From this expression 
some conclusions may be drawn: 

• If /i = At, it is useless to employ discretization 

formulae for the Riemann-Liouville derivative of 
order p higher than two because of the unavoid- 
able presence of an 0(Af)2 term. 

• The undiserable low-order term proportional to 

At is present for all WA methods with A ^ 1/2. 
Therefore, the chosen value of A (as long as A 7^ 
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Fig. 7. Numerical solution u{x, t) for the same 
problem as in Fig. 3 obtained by means of the 
weighted average method with weight factor 
A = 0.8 for 7 = 1/2, = I, S = 0.7, 
and Ax = 1/20 after 100 time steps. The 
lines are plotted as a visual aide. This case 
corresponds to the point marked by a star in 
Fig. 1. 

1/2) does not improve the precision of the WA 
method (although it matters for the stability of 
the numerical scheme). The value A = 1/2 is 
special because it removes the 0{At) term. It 
leads to the (fractional) Crank- Nicholson method. 

• The last term does not appear for non- fractional 

discretization schemes: it is characteristic of frac- 
tional methods. It becomes negligible for large 
m. In fact, for m large enough, the quantity 
^m+i^ / becomes of order of, or smaller than, 
0{At)^. The particular value of m for which this 
happens depends on the discretization formula 
of the Riemann-Liouville derivative. However, for 
the first integration steps (m small) this term, and 
consequently, the truncation error, is large unless 
the initial curvature of u{x,t), d'^u{x,t)/dx^, is 
small. These difficulties disappear for A = 1, i.e., 
they are absent in the explicit method. This sug- 
gests a practical integration procedure in which 
the first integration steps are performed by means 
of the explicit method, and the subsequent steps 
arc carried out by means of, say, the fractional 
Crank-Nicholson method. 

Of the numerical integration schemes for frac- 
tional PDE's considered here, the Crank-Nicholson 
method appears to be the most promising because 
it is always stable and is second order accurate 
in At and Ax (provided that the discretization 
formula for the Riemann-Liouville derivative is 
second order accurate in h — At). 

It has been shown that the stability of the WA 
methods can be studied by means of a von- 
Neumann-type stability analysis that is surpris- 
ingly simple and accurate. The fractional integra- 
tion schemes and stability analysis discussed here 
can be easily extended to d-dimensional fractional 
diffusive equations and fractional wave equations. 
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